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We present a modification to the Berger and Oliger adaptive mesh refinement algorithm designed 
to solve systems of coupled, non-linear, hyperbolic and elliptic partial differential equations. Such 
systems typically arise during constrained evolution of the field equations of general relativity. The 
j^jy novel aspect of this algorithm is a technique of "extrapolation and delayed solution" used to deal with 

| the non-local nature of the solution of the elliptic equations, driven by dynamical sources, within 

, the usual Berger and Oliger time-stepping framework. We show empirical results demonstrating 

the effectiveness of this technique in axisymmetric gravitational collapse simulations, and further 
t-H ' demonstrate that the solution time scales approximately linearly with problem size. We also describe 

several other details of the code, including truncation error estimation using a self-shadow hierarchy, 
■ and the refinement-boundary interpolation operators that are used to help suppress spurious high- 

frequency solution components ("noise"). 

o : 

I. INTRODUCTION 

o 

Adaptive mesh refinement (AMR) will be needed for many grid-based numerical approaches designed to solve 
a variety of problems of interest in numerical relativity, including critical gravitational collapse, binary black hole 
[ mergers, and the study of singularity structure in cosmological settings and black hole interiors [HQ, Eg- The reason 
H ' is that such problems often exhibit a wide range of relevant spatial and temporal length scales that are impossible to 
adequately resolve with a uniform mesh, given resources available on contemporary computers. In certain restricted 
bJQ' scenarios, such as the head-on collision of black holes |4}, or during the inspiral phase of a circular merger HI!, it is 
possible to construct a single, static coordinate grid that can resolve all of the length scales. However, this requires some 
a priori knowledge of the structure of the solution that will not be available in general. To date-in numerical relativity, 
mesh refinement has been used quite effectively in ID and 2D critical collapse simulations 1, 0, fToL Ull IT3 . llH , Hrf , 
the study of critical phenomena in the nonlinear sigma model in 3D Minkowski space U: , in 2D simulations of 
cosmological spacetimes , 3D simulation of gravitational waves and single black holes 17 , 0, 0, [23, 0, ^ , to 
construct initial data for binary black hole mergers [23L I24I l2o| , and the evolution of binary black hole spacetimes 

iEJ 1 . 

In the following we only consider Cauchy evolution; in other words, we have a timclikc coordinate t that foliates 
the spacetime into a set of spacelike slices, and, given initial data at t = 0, we want to evolve to the future t > 0. 
In such a coordinate basis the field equations of general relativity consist of a set of 10, second order, quasi-linear 
partial differential equations (PDEs) for 10 metric coefficients that describe the structure of spacetime. Thus, from a 
Lagrangian perspective, one would expect ten dynamical degrees of freedom per point in spacetime, where each degree 
of freedom is specified by a pair of values — a generalized coordinate and its conjugate momenta (here, for example, a 
metric element and its first time derivative). However, four of the field equations do not contain second time derivatives 
of the metric, and therefore serve as constraints, eliminating four dynamical degrees of freedom (these four equations 
are usually called the constraint equations, and the six remaining equations the evolution equations). Furthermore, 
the geometry of spacetime is invariant under arbitrary coordinate transformations of the 4 spacetime coordinates, and 
choosing a particular coordinate system (or "gauge") amounts to imposing four additional constraints, leaving only 



The notation ID refers to problems with dependence on 1 spatial dimension, in addition to the implicit dependence on time; similarly 
for 2D and 3D. 
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two dynamical degrees of freedom per spacetime point. 

The presence of constraints and coordinate freedom in the Einstein equations permits considerable leeway in the 
solution method. One of the more common methods used these days is so-called free evolution within the ADM 
( Arnowitt-Deser-Misner |28| ) decomposition (see 29] for a thorough discussion of the various possibilities and cor- 
responding classification scheme). Here, the 4 dimensional spacetime metric is written as a 3 dimensional spatial 
metric (6 independent components), lapse function and a spatial shift vector (3 components) 2 . A coordinate system 
is chosen by specifying the lapse and shift, the constraint equations are solved at the initial time, and the spatial 
metric is then evolved in time using the evolution equations. In the continuum limit, the Bianchi identities (see |27| . 
for example) guarantee that such an evolution scheme will preserve the constraints for all time, given appropriate 
boundary conditions. 

The constraint equations are elliptic in nature, whereas the evolution equations are hyperbolic. Coordinate condi- 
tions can be chosen that give algebraic, elliptic, parabolic or hyperbolic equations for the kinematical variables (lapse 
and shift). Thus, even though it is always possible in principle to adopt a free evolution approach in the numerical 
solution of Einstein's equations where elliptic equations are only solved at the initial time, there are nonetheless two 
situations where it may be preferable or necessary to solve one or more elliptic equations at each time step of the 
evolution. First, as just mentioned, it may be useful to adopt elliptic coordinate conditions. For example the choice 
of maximal slicing yields an elliptic equation for the lapse function, and the minimal distortion condition gives a set 
of elliptic equations for the shift vector components |39| . 

Second, in a numerical evolution, since the Bianchi identities will only be satisfied to within truncation error, 
the constraints can only be preserved to within truncation error 30]. This is not necessarily problematic, as the 
violation of the constraints should converge away in any consistent, stable numerical code. However, studies have 
indicated that with some formulations of the field equations, and using free evolution, certain constraint-violating 
modes grow exponentially with time, requiring prohibitively high accuracy (hence resolution) of the initial data and 
subsequent evolution to obtain a solution that is sufficiently close to the continuum one for the desired length of time 
integration^ 1^1, One possible method to circumvent this problem is to use constrained evolution rather than 
free evolution; in this case one uses some number, m, m < 4, of the constraint equations to fix m of the dynamical 
variables, in lieu of the m second-order in time equations that would otherwise be used to update those quantities 3 . 
Again, the expectation is that a stable numerical code will, in the continuum limit, provide a solution that is consistent 
with all of the field equations, including, in this case, those evolution equations that are not explicitly used in the 
overall update scheme. (Additionally, one hopes that violations in the evolution equations will not grow exponentially 
in time). To date, a significant majority of numerical codes use free evolution^], and with the exception of |38| all 
constrained evolution simulations have been carried out in ID or 2D. 4 

At this point, it is worth noting that there apparently is still an impression in the relativity community as a whole 
that elliptic equations are computationally expensive, and are thus to be avoided in numerical evolution if at all 
possible. However, provided that appropriate algorithms — such as the multigrid method — are adopted, this view 
is not consistent with at least some experience (see also 01 f° r f as t elliptic solvers using spectral methods). In 
particular, as we will show below, the solution of the elliptic equations in our adaptive code requires roughly twice 
the CPU time required to solve the hyperbolic equations; furthermore, the solution cost of the elliptics scale linearly 
with the size of the computational domain. Thus, if one considers that in a free evolution an equivalent number of 
hyperbolic equations would need to be solved in lieu of the elliptic equations, the difference in cost is not a significant 
issue in deciding whether to tackle a particular problem using constrained versus free evolution. 

We are thus lead to consider AMR algorithms for mixed elliptic/hyperbolic type, where our use of the term 
"hyperbolic" does not denote any formal definition of hyperbolicity, but instead is used to refer to an equation that 
is characterized by locality of influence. Now, a very well known AMR algorithm for hyperbolic equations is due to 
Berger and Oliger (B&O) |43_ This algorithm has several important properties that make it quite useful and efficient 
in solving certain classes of problems. These properties include dynamical regridding via local truncation error (TE) 
estimates, a grid-hierarchy composed of unigrid (single mesh) building blocks, and a recursive time-stepping algorithm 



2 We emphasize that our counting is of second-order-in-space-and-time equations. A common approach is to recast the 6 second order 
equations into a system of first order equations, which result in additional auxiliary quantities that must be evolved, and additional 
constraints among the new quantities. 

3 An alternative approach that has recently begun to receive some attention is constraint-projection (also called chopped evolution in 29]), 
where periodically during a free evolution the constraints are re-solved using a subset of t he evolved solution supplied as free data for 
the constraint solve, and afterward the evolution is continued with this new "initial data" l 371 

4 We should also mention here that there are alternative methods for evolving the field equation other than those based on a 3+1 
decomposition (see Q); in particular, character isti c or null evolution has proven useful in many situations, and recently an AMR scheme 
for characteristic evolution has been proposed l4ll . 
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that provides "optimal" efficiency in solving discretized evolution equations that are subject to a CFL-type stability 
condition. However, this algorithm, implemented verbatim for a mixed elliptic/hyperbolic system cannot be expected 
to work in general, due in part to the non-local nature of elliptic equations, and in part to the non-linear nature of 
the elliptic equations that tend to arise in numerical relativity. The reasons for this are as follows (a more detailed 
discussion is given in Sec. IHII and Sec. II V|) . In the B&O time-stepping procedure, a single, large time step is taken on 
a coarser level before several smaller time steps are taken on the interior, fine level. This is done so that the solution 
obtained on the coarse level can be used to set boundary conditions (via interpolation in time) for the subsequent 
fine level evolution. As the hierarchy is generated via local truncation error estimates, the solution obtained on the 
coarse level in the vicinity of the fine level boundaries will presumably be sufficiently accurate to allow one to use 
the coarse level solution to set fine level boundary conditions without adversely affecting the global solution. If the 
equations are hyperbolic, then a poorly resolved solution in the interior region of the coarse level will not have time 
to "pollute" the coarse/fine boundary region in only a single coarse level evolution step (and the coarse level solution 
is refreshed every time step with the fine level solution in the injection phase of the algorithm). This last statement 
is not true for elliptic equations in general, for then poorly resolved source functions in the interior of the coarse level 
could globally affect the accuracy of the solution obtained during the coarse level evolution step. For certain kinds 
of linear elliptic equations, such as the Poisson equation in Newtonian gravity, or that arising in the incompressible 
Navier-Stokes equations, one can circumvent this problem by taking advantage of conservation laws satisfied by source 
fields that couple to the elliptic equtions (for example, with the Poisson equation in Newtonian gravity one can use a 
fine to coarse level injection function that preserves the matter energy density, and see for example [31| for a method 
for solving the Navier-Stokes equations with B&O style time sub-cycling). In general relativity the equations are 
non-linear, and furthermore are coupled in such a manner that it is impossible to isolate such source functions in 
general. Therefore, to use B&O AMR in a constrained evolution, in particular its time-stepping algorithm, requires 
some modifications to deal with the elliptic equations; these modifications are the prime focus of this paper. 

In Sec. m we first review the axisymmetric gravitational collapse code introduced in |32j that was used to develop 
the AMR approach described in this paper. This code solves a discretized version of the Einstein-Klein Gordon system 
of equations. We then review the original Berger and Oliger algorithm in Sec. IIIII and then proceed to a description 
of the modifications we have made to handle elliptic equations in Sec. IIVI Our primary modification involves a split 
of the solution of the elliptic equations into two phases. During the first phase, when hyperbolic equations are solved, 
functions satisfying elliptic equations are extrapolated to the advanced time level. The second phase is delayed until 
all finer levels have been evolved in the same fashion via recursion, and hence all levels at the same or finer resolution 
as the given one are in sync (i.e. have been evolved to the same physical time). Then, the elliptic equations are 
solved over the entire sub- hierarchy from the given level to the finest, using extrapolated boundary conditions from 
the parent (coarser) level at interior coarse-grid boundaries. In Sec. we present several simulation results, including 
convergence tests and comparison with unigrid simulations. Concluding remarks are given in Sec. I VII All equations 
and finite-difference operators are listed in App. ^ and some additional details of the AMR algorithm are given in 

App. m 

II. AN AXISYMMETRIC GRAVITATIONAL COLLAPSE CODE 

In this section we briefly review the physical system we are modeling (general relativity with a scalar field matter 
source), the PDEs governing the model, and the unigrid numerical code that computes an approximate finite-difference 
solution of the PDEs; additional details can be found in [33. |43| . 

A. Equations, Coordinate System and Variables 

The Einstein field equations can be written as 

Rnu - -^Rg^u = SnT^, (1) 

where is the Ricci tensor, R = is the Ricci scalar (using the Einstein summation convention where repeated 
indices are summed over), and we use geometric units where Newton's constant, G, and the speed of light, c, are set 
to 1 . With a massless scalar field $ as the matter source, the stress-energy tensor Ta V is given by 



(2) 
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and the evolution of <E> is governed by the wave equation 

□$ = $./ = 0. (3) 

In these expressions a comma (,) is used to denote a partial derivative, and a semicolon (;) a covariant derivative. 

Restricting attention to axisymmetric spacetimes without angular momentum, and choosing cylindrical coordinates, 
(t, p, z, <p), adapted to the symmetry, we can write the spacetime metric as 

ds 2 = -~a 2 dt 2 + ?/> 4 [ (dp + (3 p dt) 2 + (dz + (3 z dtf + p 2 e 2pS d<j) 2 ] . (4) 

The axial Killing vector is (d/dcj))^ and hence all the metric functions a, j3 p , (3 Z ,ip and er, and the scalar field $ depend 
only on p, z and t. Almost all coordinate freedom has been eliminated by choosing this form for the metric. What 
remains to be specified is a time-slicing, and for this we use maximal slicing, defined by 

dK 

K = 0, = 0, (5) 

where K = K a a is the trace of the extrinsic curvature tensor K a b of t — const, slices 39]. This gives an elliptic 
equation for a (|A2(I . The constraint equation subset of Q gives 3 additional elliptic equations: the Hamiltonian 
constraint, which is viewed as an equation for -0 (|A3(I . and the p and z components of the momentum constraint, 
which are treated as equations for f} p (|A4|) and /3 Z (|A5ll , respectively. One member of the evolution subset of yields 
a second order evolution equation for a, and the wave equation (jSJ) provides a second order hyperbolic equation for 
$. We convert both of these evolution equations to first-order-in time-form l|A7IA9(l by defining conjugate variables 
f2 (geometry) and II (matter) as follows: 



-a. t + 2{3 p (pa) , + f3 z a, : 



P 



(G) 



n = — ($ t -/3 p $.p + /3 2 $, z ). 

a 



(7) 



Thus we end up with a system of 8 equations for 8 variables — a,ip,(3 p and (3 Z satisfy elliptic equations, and a, fi, i> 
and II satisfy hyperbolic equations. 



B. Boundary Conditions 

On the axis at p = 0, the following regularity conditions must be enforced in order that spacetime remain locally 
flat in the vicinity of the axis: 

a P = 0, 
t P = o, 
/3,p = 0, 

pp = 0, 

(7 = 0, 

= 0, 
<f.p = o, 

n P = o. (8) 

At the outer boundaries p — p ma x, z — z max and z = z m i n , for the hyperbolic variables approximate outgoing 
radiation (Sommerfeld) conditions are imposed 5 , while for the elliptic equations, conditions based on asymptotic 



5 These conditions assume that spacetime is nearly flat at the outer boundary, and that disturbances (radiation) in both the scalar and 
gravitational fields are propagating purely radially, and have 1/r fallofT. 



5 



As an initial guess to the solution at time ti , copy variables 
from to to t\ ; 



repeat 

perform 1 Newton-Gauss-Seidel relaxation sweep of the 

evolution equations, solving for the unknowns at time t\ ; 

perform 1 multigrid vcycle on the set of elliptic equations, 
discretized at time t\ ; 
until (residual norm < tolerance) 



FIG. 1: A pseudo-code description of the iteration we use to solve the system of coupled hyperbolic/elliptic equations on a 
single mesh. 



flatness conditions are used: 



a = 1, 

ip — l + p^.p + zi>,z = o, 

(3 Z = 0, 

d" = 0, 

- zo z + a = 0, 

ztt tZ + = 0, 

z$ 2 + $ = 0, 

rn t + pn „ + z n z + n = o. 



r$ t - 



(9) 



C. Unigrid Numerical Scheme 



The set of 8 PDEs are solved using second-order accurate finite difference (FD) techniques. The elliptic equations 
are solved using an FAS (Full Approximation Storage) multigrid algorithm with ^-cycling 0> E5- At t = 0, 
ct, and II are freely specified, after which the remaining variables a,(3 p ,/3 z and ip are obtained by solving the 
corresponding elliptic equations. To evolve the variables with time, the hyperbolic equations are discretized using 
an iterative Crank- Nicholson scheme, with Kreiss-Oliger 44J dissipation added to reduce unwanted (and un-physical) 
high-frequency solution components ("noise") of the FD equations. This iteration involves variables at two time levels: 
the known solution at t — t$, and the unknowns, solved for using Newton-Gauss-Seidel relaxation implemented in 
RNPL |45(, at t = t\ = to + At. After each iteration, the elliptic variables are updated at t = t% by applying a 
single I^-cycle. This process is repeated until the infinity norm of the residual of all equations is below some specified 
tolerance. The pseudo-code in Fig. ^ summarizes this iteration sequence: 

Specific difference operators used to discretize the equations are summarized in App. IA 31 



III. THE BERGER AND OLIGER AMR ALGORITHM 



Here we briefly review some aspects of the B&O AMR algorithm for hyperbolic PDEs that are of relevance to this 
paper, in particular the grid hierarchy and time-stepping procedure. 



A. AMR Grid Hierarchy 



In the Berger and Oliger AMR algorithm, the computational domain is decomposed into a hierarchy of uniform 
meshes (see Fig. |2J) with the following properties: 

• The hierarchy contains If levels. Each level £ contains grids of the same resolution — the coarsest grids are in 
level 1 (the base grid), the next-coarsest in level 2, and so on until level £/, which contains the finest grids in 
the hierarchy. 
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• The ratio of discretization scales hgjhg+x between levels £ and £ + 1 is called the spatial refinement ratio p s ^. 
p St i is typically an integer greater than or equal to 2. For simplicity, we will also assume that p s i is the same 
for all levels, and therefore use the symbol p s to denote the spatial refinement ratio. 

• All grids at level I + 1 (child grids) are entirely contained within grids at level £ (parent grids). Grids at the 
same level may overlap. 

• In the simplified variant of the B&O algorithm described here, we require that all grids within the hierarchy share 
the same coordinate system. In particular, this implies that all grid boundaries run parallel to the corresponding 
boundaries of the computational domain. In addition, a child grid must be aligned relative to its parent grid 
such that all points on the parent grid, within the common overlap region, are coincident with a point on the 
child level. The original B&O algorithm allowed for a child grid to be rotated relative to its parent. 

The particular grid structure that exists at any given time is calculated by computing local truncation error 
estimates, so that at any point x = (p, z) within the computational domain the finest grid covering that point has 
sufficient resolution to adequately resolve all features of the solution there. This is an important property of the grid 
hierarchy, not only for the obvious reason of providing the desired resolution everywhere, but it justifies the use of 
the B&O time-stepping algorithm to evolve the hierarchy, as we now review. 

B. The Berger and Oliger Time-Stepping Algorithm 

The B&O recursive time-stepping algorithm was designed to solve hyperbolic equations, discretized on the AMR 
grid hierarchy. The basic ideas behind this update scheme are as follows. The hierarchy is evolved in time through a 
particular sequence of unigrid time-steps, performed on individual grids within the hierarchy. A time step of size Atg 
is taken on all grids at level £, before a number p t ^ (the temporal refinement ratio) time steps of size Ai^ + i = Ate/pt,i 
are taken on level i + 1. The preceding rule is applied recursively, from the coarsest to finest level in the hierarchy. 
In general, pt/ must be an integer greater than or equal to p s ^ in order to satisfy the CFL condition on all levels in 
the hierarchy if it is satisfied on the coarsest level. As with p s , we only consider a constant temporal refinement ratio 
Pt for all levels. The reason why a time step is first taken on coarse level £ is that the solution obtained there at time 
t + Ate is then used to set boundary conditions, via interpolation in time, for the subsequent time steps on the finer 
level £ + 1 (unless some portion of the fine level abuts the boundary of the computational domain, in which case the 
physical/mathematical boundary conditions of the original problem can be applied.) It is possible to do this because 
the solution obtained on the coarse level in the vicinity of the finer level boundary will be as accurate, to within the 
specified truncation error, as a putative solution would have been that was obtained on a fine level encompassing 
the entire computational domain 6 . The solution obtained on the coarse level interior to this boundary will not be 
as accurate; however, the assumed hyperbolic nature of the PDEs will protect this inaccuracy from polluting the 
coarse/fine boundary region within a single coarse level time step. 

After p t time-steps on level £ + 1, when the solution on grids at levels £ and £ + 1 are again in synchrony, grid 
functions from level £ + 1 are injected into the coarse grids at level £, in the region of overlap between the two levels. 
Thus, the most accurate solution available at a given point x is continuously propagated to all grids in the hierarchy 
that contain x. Injection simply consists of copying values from level I + 1 to level I at common points (in the more 
general B&O algorithm, where finer levels can be rotated relative to coarser levels, the injection step requires some 
form of interpolation). Fig- El is a pseudo code description of the B&O time-stepping procedure just described. 

IV. BERGER AND OLIGER STYLE AMR FOR CONSTRAINED EVOLUTION 

The locality argument given in the preceding section to justify the use of the B&O evolution algorithm is only 
applicable to hyperbolic equations, for then the finite speed of propagation prevents contamination of the solution in 
the boundary region of a coarse level by a poorly resolved solution in the interior. Solutions to elliptic equations do 
not share this property, and therefore it is not feasible to solve for such equations on the coarse grid alone, with the 



Assuming that the specified maximum TE estimate used in the construction of the grid hierarchy is sufficiently small that the solution 
is within the convergent regime, and hence the TE estimate is a good approximation to the actual solution error. 
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Computational domain, 
covered by a hierarchy 
of uniform grids 



Level 1 grid 



Level 2 grids Level 3 grids 




FIG. 2: An example of a Berger and Oliger mesh hierarchy. The hierarchy consists of a number of levels, where each level 
contains a set of uniform meshes of the same spatial resolution. Our convention is to let higher level numbers denote levels 
consisting of grids with higher resolution (smaller mesh spacing, h), beginning with level 1 for the coarsest mesh. The upper 
diagram shows the computational domain, covered by a three-level-deep hierarchy. The plots below this demonstrate how the 
hierarchy is stored in memory, namely as a collection of individual grids. Thus a given point, x in the computational domain 
can be contained/represented in multiple grids in the Berger and Oliger scheme. 

intention of supplying boundary conditions for subsequent fine grid time steps 7 . One way to circumvent this problem 



7 Except, as mentioned in the introduction, for certain kinds of linear systems, where the source terms appearing in the elliptic equations 
can cleanly be identified and coarsened in a manner that conserves the source "energy density". This is not possible for the Einstein 
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repeat N times: call single_step(l) ; 
stop; 

subroutine single_step (level £) 

if (regridding time for level £) then 

regrid from levels i to if; 
end if 

if {£ >1) then 

set boundary conditions along AMR boundaries at time 
tt + Ati via interpolation from level (£—1); 

end if 

perform 1 evolution step on all grids at level £; 

te, ■= tg + Atg ; 

if (£ < £ f ) then 

repeat [p t := Aif/Ai^+i] times: call single_step(£ + 1) ; 
inject the solution from level i+l to level t 

in the region of overlap between levels £ and £ + 1 ; 

end if 

end of subroutine single_step 

FIG. 3: A pseudo-code representation of the Berger and Oliger time stepping algorithm. 

is to abandon the B&O recursive time-stepping procedure. In other words, one could evolve the entire hierarchy 
forward in time with a global time step, for example by performing a Crank- Nicholson style iteration as in the unigrid 
code (see Sec. Ill C|) . A major drawback to this method is that, to satisfy the CFL condition, the global time step will 

need to be set to \Ap(£f), where Ap(£f) is the cell size on the finest level £f in the hierarchy, and < A < 1 is a 

it— i 

constant. This would require that one take p t — 1 additional time steps at level £ for each time step that the usual 
BhO algorithm would have taken. 

The technique that we propose here to incorporate elliptic equations into the standard B&O time-stepping frame- 
work is to employ a combination of extrapolation and delayed solution of the elliptic variables (see Fig. [3J. Simply 
stated, on coarse levels one does not solve the elliptic equations during the evolution step of the algorithm; rather, 
one extrapolates the corresponding variables to the advanced time from the solution obtained at earlier times. The 
solution of the elliptic equations is delayed until after the injection of fine grid (level £ + 1) values into the parental 
coarse grids (level £). At that stage, all levels from £ to £/ are in sync, and the elliptic equations are solved over the 
entire resulting subset of the hierarchy, with boundary conditions on AMR boundaries of level £ set via extrapolation. 
This ensures that all details from finer grids interior to level £ are represented in the solution. 

One of the non-trivial aspects of this technique is the method used to extrapolate the elliptic variables, which we 
now describe. We use linear (2nd order) extrapolation in time, with periodic corrections to try to account for changes 
that occur upon global multigrid solves. For level £, £ > 1, the two past-times used in the extrapolation are the 
two most recent times when levels £ and £ — 1 were in sync (thus, at times when a solution of the elliptic variables 
involving at least levels £—1 to £f was obtained); in other words, every p t steps we save the elliptic variables for use in 
extrapolation 8 (see Sec. IB 51 for a pseudo-code description of how the past time levels are initialized for the very first 
time step of evolution). For level 6 = 1, the two most recent time levels are used for extrapolation 9 . The correction, 
applied only to levels £ > 1, is calculated as follows. Whenever level £ is in sync with level £ c ,£ c + l, ...,£/ + !, £f, where 



equations, which are fundamentally non-linear, since the gravitational field acts as its own source. 

Early experiments using the values from the two most recent time steps of level I for extrapolation resulted in unstable evolution. 
Though note that in our algorithm level 1 is always fully refined because of the self-shadow hierarchy mechanism we use for truncation 
error estimation (see Sec. IB II . and therefore level 2 should be considered the "true" coarsest level of the hierarchy (in effect, level 1 is 
only used for truncation error estimation) 
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repeat N times: call single_step(l) ; 
stop; 

subroutine single_step (level £) 

if (regridding time for level £) then 

regrid from levels i to if; 
end if 

for hyperbolic variables: if (£ > 1) then set boundary conditions along 
AMR boundaries at time (tg + Atg) via interpolation from level £—1; 

for elliptic variables: extrapolate the entire grid function to time 
(ti + Ati) from earlier-time values; 

repeat 

if (£ — £f) then perform 1 multigrid vcycle on elliptic variables 

at time (ti + Atg) ; 
perform 1 iteration of the Crank-Nicholson sweep for all 
hyperbolic variables; 
until (residual norm < tolerance) 

ti := ti + At£ ; 

if (£ < £ f ) then 

repeat [pt := Ai^/At^+i] times: call single_step(£ + 1) ; 
compute the truncation error estimate for level £ + 1 by subtracting 
the solution in the region of overlap between levels £ and £ + 1 ; 
inject the solution from level £+1 to level £ 

in the region of overlap between levels £ and £ + 1 ; 

end if 

if (£=1 or {£<£ } and t(\=t(,-i)) then 
t := ti; 

do £ -.= 1 + 1 to £ f 

for each elliptic variable ft (t) : fi (t) ■= ft {t); 
end do 

re-solve the elliptic equations over the sub-hierarchy at time t; 

do £ := £+1 to £ f 

for each elliptic variable fe (t) : 

A/^W :=/*,(*)- 4 (*); 

/ c , (i) :=A/, (t)/^; 

// (t - p t At, ) := fi (t - ptAtto) + Af e „(t) - f cto (t); 

end do 
end if 

end of subroutine single_step 



FIG. 4: A pseudo-code representation of the modified Berger and Oliger time stepping algorithm described in Sec, II VI compare 
to the original algorithm in Fig. [3] Notice that here we have expanded the "perform 1 evolution step..." statement in Fig. 
to highlight the fact that in the modified algorithm the finest level is treated differently than the coarser levels then (though 
the particular details of the evolution step are not significant — for concreteness we list the same scheme as used in the unigrid 
code). Here we also show where the truncation error estimate is computed when using a self-shadow hierarchy (see App. 15} . 
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£ c < £, a multigrid solve takes place over levels £ c ,...,£f. Denote by ft(t) the value of a variable / at level £, calculated 
via extrapolation from f(t — p t Ati)i and f(t — 2p t Ati)i 1 and let fg(t) denote the value of the same variable after the 
multigrid solve at time t (note that for simplicity in notation we have dropped the spatial coordinate dependence of 
the variable /). As illustrated in Fig[3 the correction contains two components: 

Aft{t) = f e (t) - h{t) and / c( (()e^P, (10) 

Pt c 

which are used to change the past time value ft (t — p t Ati) as follows: 

f e {t - Pt At e ) — > f £ (t - p t Ati) + Afe(t) - f ct (t) (11) 

The logic behind this form of correction stems from a couple of observations about how the re-solved solution differs 
from the extrapolated solution, and how some adaptive solutions differ from unigrid solutions of comparable resolution. 
First, in general Afi(t) is (in a loose sense) proportional to £ c — I; i.e. the more levels over which the elliptic equations 
are re-solved, the larger the change in the interior, finer level t solution. However, the change in the interior part 
of the solution induced by the global solve tends to be a near constant shift, leaving finer details of the interior 
solution unchanged. Second, the local "velocity" fi{t) — f(t — p t Ati)i (calculated prior to the correction) tends to be 
represented quite accurately in the adaptive solution scheme. Therefore the Aft(t) part of the correction preserves this 
velocity for subsequent extrapolation, and is in fact essential for the stability of the algorithm with a deep hierarchy, 
as then the global shift is often larger in magnitude than the local velocity. The second part of the correction f c e(t) is 
an attempt to improve the accuracy of the extrapolation, for if in hindsight the quantity f c i(t) had been added to / 
at each one of the p\ ~ lc intermediate time steps between solves over levels £ c --£f, then Afi(t) would be zero at time 
t (ignoring, of course, the effect that this putative correction would have had on the solution). 

One last comment regarding the extrapolation: on the finest level of the hierarchy, the elliptic equations are solved 
via the usual interleaved Crank- Nicholson/F-cycle iteration as discussed in Sec. Ill Cl then the extrapolation simply 
serves to set boundary conditions. 



V. RESULTS 



The extrapolation technique described in the previous section is rather ad-hoc, though both the choice of which 
past times to extrapolate from and the use of corrections play a significant role in the stability and accuracy of the 
adaptive evolution code. In this section we present some results showing the effectiveness of the algorithm. 



A. Comparison with a Unigrid Evolution 



The first test of the algorithm presented here is a comparison of unigrid evolutions to a similar AMR evolution 
(test results on the convergence properties of the unigrid code can be found in |32|'). Specifically, we compare an 
evolution obtained with the AMR code to a unigrid evolution, where the entire unigrid mesh is given the resolution h 
of the finest level in the AMR hierarchy. To gauge how well the AMR solution approximates the unigrid one, we then 
compare the unigrid solution to that obtained with two additional unigrid runs, with resolutions of 2h and Ah. In 
certain respects this is not a very stringent test, as limited computational resources do not allow us to run using very 
high resolution unigrid simulations (which of course is the motivation for pursuing AMR). Also, the accuracy of the 
AMR solution will largely depend on the structure of the grid hierarchy, which is (predominantly) controlled by the 
maximum allowed TE. Thus, in principle one could obtain better agreement between the AMR and unigrid simulations 
by decreasing the TE parameter while keeping the maximum depth of the AMR hierarchy fixed. Nevertheless, this 
comparison does demonstrate that the adaptive algorithm works in practice. 

The initial data for this example is a time symmetric scalar field pulse: 



$(0, p, z) — Aexp 



p 2 + z 2 



(12) 



with A — 0.23 and A = 1.0 (all other free fields are set to at t — 0). This amplitude is sufficiently large 
that the solution is in the non- linear regime (and thus close to forming a black hole). The outer boundary is at 
Pmax = z max = — z min = 10. The maximum value for the TE is set to 10 -5 (only $ and II are used in the calculation 
of the TE — see App. [Bj. The resolution of the base grid is 65 x 129, with p s = 2, and up to 3 additional levels of 
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t - AAt t - 3At t - 2At t- At t 




t - 4At t - 3At t - 2At t - At t 




t - 4At t - 3At t - 2At t- At t t + At 



FIG. 5: An illustration of the technique used to extrapolate and solve for elliptic variables within the AMR framework. In this 
example, we show the evolution of an elliptic variable fa from t — At to t, assuming that p t = 2 and I < if. Initially (step 1), 
fa at t is calculated via linear extrapolation from data at past times t — 2At and t — 4At (quantities used for extrapolation 
are depicted by boxes in the diagram). This value, labeled fait), is used during the solution of the hyperbolic variables (via 
Crank-Nicholson iteration (CN) here). At time t, we assume levels £ and £ c = t — 1 are in sync, and so after the CN iteration 
(and after the equations on finer levels £ + l..£f have been evolved to t) the elliptic equations are re-solved over levels £ c ..£f 
(step 2). This results in a change in the value of / from fa{t) to fe(t) (for simplicity we assume that a similar change that 
occurred at time t — 2At is zero). This change is propagated back to t — 2At, so that the same velocity v, modulo a correction 
fci{t), will be used to extrapolate / to t + At (step 3). Here, since £ — £ c = 1, the correction is such that one is effectively 
extrapolating from t — A At and t to time t + At; this would not be the case otherwise — see (I1UI1H . (If £ = if , then the elliptic 
variables are solved within the CN iteration, and the value obtained afterward is used for fa(t).) 

refinement. Thus, level 4 has an effective resolution of 513 x 1025, and so we choose unigrid comparison runs with 
resolutions of 513 x 1025 (h), 257 x 513 (2h) and 129 x 257 (4h). The Courant factor is A = 0.3 for all runs, and the 
physical time at the end of each simulation is t ~ 2.8 (which corresponds to 480 time steps on the level with the finest 
resolution) . 

Fig. shows a plot of the conformal factor rp at t = 2 from the adaptive simulation, with grid bounding boxes 
overlaid to give an idea of the structure of the hierarchy adaptively generated by the AMR code. Fig. [7\ shows £2 
norms of the differences between the solutions generated by the h resolution unigrid simulation and the two lower 



12 




FIG. 6: A surface plot of ijj at t = 2 from the adaptive simulation discussed in Sec. IV Al where the height of the surface is 
proportional to the magnitude of tj) (ranging approximately from 1 at the outer boundary to 1.5 at the origin p = z = 0). 
Overlaid on the surface are the AMR grid bounding boxes — the smallest, interior box has the finest resolution, corresponding 
to an effective unigrid resolution of 513 x 1025. 

resolution unigrid and adaptive simulations. For brevity we only show differences for the 4 elliptic variables; differences 
in the hyperbolic variables exhibit similar behavior. What Fig. \7\ demonstrates is that, compared to the h unigrid 
simulation, the adaptive solution has accuracy comparable to the 2h simulation, and significantly greater accuracy 
than the 4/i simulation, even though (as illustrated in Fig. [|JJ| the majority of the coordinate domain in the adaptive 
solution is spanned by a grid with less resolution that either that of the 2h or Ah simulations. Fig. [S] is a comparison 
of the time-difference of the maximum (or minimum, as appropriate) of the elliptic variables ij), a, (3 P and f3 z from the 
h unigrid and adaptive simulations. The time-difference A// At is calculated as (f n+1 — f n )/At, where /" denotes 
one of the quantities, ip,a,f3 p or (3 Z at time step n. Fig. [S] demonstrates two interesting aspects of the extrapolation 
scheme of the adaptive code. First, the high frequency temporal "noise" that is apparent in the adaptive solution of 
the elliptic variables (and that has been exaggerated in the figure by taking a time-difference) does not adversely affect 
the accuracy of the solution on average. Second, the presence of such noise suggests why linear extrapolation at level 
I from the most recent times of level i is unstable; however, it is not so obvious why extrapolation from past times 
that are in-sync with a parent level (£ — 1 or less) results in stable evolution. Note that with the particular system 
of equations solved for here (summarized in App. [SJ no explicit time derivatives of any elliptic quantities appear. 
Furthermore, during the evolution phase of the algorithm (see Sec. Ill C|) . the Crank- Nicholson differencing scheme only 
couples the time average (over two time steps) of elliptics variables to the hyperbolic equations, providing a certain 
amount of temporal smoothing. Thus, the results shown in Fig. [S] suggest that modifications to the extrapolation 
scheme may be needed for systems of PDEs that use other methods to difference in time, or if time derivatives of 
variables solved for using elliptic equations couple to the hyperbolic equations. 

1. Timing information 

We conclude this section by presenting timing information in Tab. [I] below for the comparison test just described. 
This serves to show that, at least for the elliptic- hyperbolic system considered here, the overhead of the AMR algorithm 
is negligible, and that the solution time scales roughly linearly with the total number of grid points in the discrete 
solution. See also Sec. IV B II which contains more detailed information on the percentage of time spent in various 
stages of the algorithm for the test simulations presented in the next section. 
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FIG. 7: Comparison between unigrid and adaptive simulation results for the elliptic variables if),a,l3 p and /3 Z , as discussed 
in Sec. IV Al Shown here are £2 norms of differences between the solution generated by the h resolution (513 x 1025) unigrid 
simulation and each of the solutions produced by two lower resolution unigrid simulations, 2h (257 x 513) and 4/i (129 x 257), 
and an adaptive simulation (see Fig.|S]for a representative sample of the mesh structure from the AMR solution at t = 2, where 
the base level 1 has resolution 65 x 129, and the finest level 4 has the same resolution as the h unigrid run). Thus, we are 
using the h unigrid solution as a benchmark, and the plots show that the adaptive solution is of comparable accuracy to the 2h 
unigrid solution, and of significantly greater accuracy than the 4/i unigrid solution, despite that the majority of the coordinate 
domain of the AMR solution is covered by grids with mesh spacing hi satisfying hi > 4/i > 2h. 



B. Convergence and Consistency Tests 



In this section we present some convergence test results, providing evidence that the overall AMR solution scheme 
is stable and convergent. We also compare the solution of the conformal factor xj} to that obtained in a partially 
constrained evolution, where a,/3 p and (3 Z are solved as described in Sec. [H] via the slicing condition and momentum 
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FIG. 8: Comparison between h resolution (513 x 1025) unigrid and adaptive simulations results as discussed in Sec. IV Al 
showing the time-difference of the maximum (or minimum, as appropriate) of the elliptic variables ip, a, f5 p and /3 Z . For clarity, 
only about 50 time-steps are shown, and the differences for a and ip have been scaled by constants to fit all the plots on the 
same vertical scale. The time-difference Af/At is calculated as (f" +1 — f n )/At, where f n labels one of the above quantities 
at time step n. These figures demonstrate a couple of interesting aspects of the extrapolation scheme for the adaptive code. 
First, the high frequency "noise" that is apparent in the adaptive solution (which has been exaggerated here by taking a 
time-difference) does not adversely affect the accuracy of the solution on average. Second, the presence of such noise suggests 
why linear extrapolation from the most recent time levels is unstable (though it is not obvious why extrapolation from past 
time levels that are in-sync with a parent level results in stable evolution). 

constraints respectively, but where a hyperbolic evolution equation (|A10|I obtained from the maximal slicing condition 
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umgrid-4n 


unigrid-zn 


unigrid-n 


AMR 


runtime (s) 
total number of grid points 
average time per grid point (fis) 


4.45 x 10 2 
4.01 x 10 6 
111 


4.48 x 10 3 
3.18 x 10 7 
141 


3.31 x 10 4 
2.54 x 10 8 
131 


2.19 x 10 3 
1.97 x 10 7 
111 



TABLE I: Timing information for the tests described in Sec IV Al The runtime is the wall time taken (on a 2.4Ghz Intel XEON 
processor), in seconds, for each simulation, including initial data calculation and evolution. The total number of grid points is a 
count of all the grid points, in space and time, at which a solution was obtained during the simulation. This includes the grid 
points used during calculation of the initial data, and the calculation of the initial hierarchy for the adaptive run. The average 
time per grid point (in microseconds) is the runtime divided by the total number of grid points. What these numbers suggest 
(see also Tab. [n]is that the computational cost of this solution method scales approximately linearly with the total number of 
grid points, and that the computational overhead for the adaptive algorithm is negligible. Note however that we have not taken 
into account that in the adaptive hierarchy certain coordinate locations are covered by multiple points, and therefore there is 
some "wastage" in the AMR calculation. In this particular simulation, at any one time between 15 — 20% of the points in the 
adaptive hierarchy were redundant, and this should be viewed as an additional (and unavoidable) overhead to the per-point 
cost of this Berger and Oliger style AMR algorithm. 



K = 0, rather that the elliptic Hamiltonian constraint, is now used to update ip 10 • in t ne continuum limit, and the 
limit where the outer boundary position goes to infinity, the solution obtained from fully and partially constrained 
evolution should (assuming both numerical implementations are consistent and stable) converge to a unique solution 
of the Einstein equations. In the partially constrained evolution, ip is evolved as a hyperbolic variable, and thus within 
the traditional B&O time stepping framework; however, during fully constrained evolution ip is solved for using an 
elliptic equation using the new AMR technique. Thus, demonstrating convergence to a consistent solution for ip is a 
rather non-trivial test of the modified AMR algorithm. 

We use the following technique to calculate convergence factors for the adaptive code. We choose a "modest" value 
(10 -3 in this case) for the maximum allowed TE (calculated as described in Sec. IB 1|) . run a simulation, and save a 
copy of the dynamical grid structure produced during the evolution. The solution on this grid hierarchy will serve as 
the coarsest resolution simulation, labeled Ah. For the higher resolution simulations 2h and h, We rerun the code with 
identical initial data, and use the same grid hierarchy produced by the Ah case except that the resolution of all grids 
in the hierarchy is doubled (quadrupled) for the 2h (h) simulation; the Courant factor is kept constant, hence the 
number of time-steps taken per level also doubles (quadruples). Then, by assuming the usual Richardson expansion, 



one can compute a convergence factor Q h f for a variable / via 



nh _ \\fah - hh\\2 , 1Q s 

\\j2h - fh\\2 

where fih is the solution on the Ah hierarchy, and similarly for and //,.. In 113|) . the subtraction of grid functions 
is only defined at common points between the hierarchies (every other point of the finer mesh in this case) . The £2 
norm of a grid function / is taken point-wise as follows 

1/2 

II = v - 9 l - J '- — . ( I -I ) 

EEE 1 

l g i,j 

where the sum over I is the sum over all levels in the hierarchy where / is defined (so for the differences in (|13|l this 
will be all levels except the finest), the sum over g is over all grids at the given level, and the sum over (i,j) covers all 
points within a given grid, where {p g , z g ) is the location of the (i = 0, j = 0) point. Note that such a point-wise norm 
is biased (compared to an area-weighted norm, such as an integral norm) toward the highest resolution region of the 
domain, where most of the points are clustered. This is desirable for the solution presented here, where the region 
of high refinement is centered on a very small part of the domain, and outside of this region the solution is slowly 



That an independent equation (other than the Hamiltonian constraint) exists for tp is due to the over-determinism in the Einstein 
equations, as discussed in the introduction. 
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AMR-4h 


AMR-2h 


AMR-h 


runtime (s) 


3.38 x 10 2 


2.00 x 10 3 


1.40 x 10 4 


percentage of runtime solving elliptics equations 
percentage of runtime solving hyperbolic equations 
percentage of runtime in AMR related functions 
percentage of runtime in miscellaneous functions 


68.2% 
24.4% 
3.4% 
4.0% 


60.3% 
33.1% 
3.4% 
3.2% 


50.8% 
43.1% 
3.8% 
2.3% 



TABLE II: Timing information for the fully constrained simulations described in Sec lVBI The runtime is the wall time taken 
(on a 2.4Ghz Intel XEON processor), in seconds, for each simulation, including initial data calculation and evolution. For this 
set of simulations, the total number of grid points in space and time (compare Tab. HJ increases by a factor of 8 from the 4/i to 
2h simulation, and again by a factor of 8 from the 2h to h simulation. Hence the runtime should increase by a corresponding 
factor for an algorithm whose cost scales linearly with the number of points — we do, approximately, see this behavior. The next 
four rows give a rough breakdown of the percentage of total time spent in each of the key parts of the program (calculated using 
the function profiling option of Portland Group's pg compilers). The two main points we want to stress with this data is that 
the scaling (with problem size) of the adaptive multigrid algorithm used to solve the elliptic equations is close to linear, and 
is not significantly slower than the solution of the hyperbolic equations. Second, the cost of AMR related functions (including 
regridding, truncation error estimate calculation, interpolation and injection functions) is quite small compared to the cost of 
solving the numerical equations. The miscellaneous functions of the last row include calculations of diagnostic quantities, the 
cost of the apparent horizon finder (that searches for the presence of black holes), and I/O. 



varying and well represented by the coarse mesh. An area-weighted norm in this case would almost completely ignore 
information on the finest levels. For a "residual" function R, in other words a function that should converge to zero 
in the limit as h — > 0, we use the following expression for its convergence factor: 

nh _ \\R2hW2 n 

°* = P^F (15) 

For a second-order accurate finite difference scheme one expects both Q*j and Q\ to approach 4 asymptotically. 

The initial data for this example is also a time symmetric scalar field pulse given by (|12fl . with A — 0.25 and 
A = 0.5, and all other free fields set to at t = Q. Adopting apparent horizon detection as our operative definition of 
black hole existence, in this case a black hole (with mass M ~ 0.12 in geometric units) does form by t ~ 2.5. 

Computationally, relevant parameter settings are as follows. The outer boundary is at p max = z max = — Zmin = 32, 
and for the Ah simulation the resolution of the base grid is 65 x 129. The maximum value for the TE is set to 10 -3 ; 
this results in a grid hierarchy containing 3 additional levels (p 8 = 2 : 1) at t = 0, and 6 additional levels at the end of 
the simulation 11 . See Fig.[5]for sample plots of 4> at t = and t = 3 to illustrate the grid hierarchy. The effective finest 
grid resolution for the h simulation is roughly 16, 000 x 32, 000, making it impractical to do a unigrid comparison in 
this case. Fig. 1101 shows the calculated convergence factors for the four elliptic quantities ip, a, (5 P and f3 z , and Fig. 1111 
contains the convergence factor for ip c — tpf, where ip c is the conformal factor ip from fully constrained evolution, 
and ipf is the conformal factor calculated from the free evolution of if). These plots show reasonable convergence and 
consistency results, with a couple of caveats discussed in the captions. 



1. Timing information 

Tab. [H] contains some timing information for the set of simulations described in this section. The major point we 
would like to emphasize with this table (see the caption for more information) is that solution of the elliptic equations 
is not significantly more expensive than the solution of the hyperbolic equations in our model; moreover the scaling of 
the solution time with problem size is close to linear, again for both types of equations. This is not surprising or new 
in any sense, as one can immediately predict this "Golden Rule" scaling behavior from Brandt's work on multi- level 
adaptive (MLAT) schemes 0|; indeed, that elliptic equations could be solved in linear time within the context of 
general relativity was already demonstrated in the early 1980's 0,0]. 



11 the increase in hierarchy depth that occurs when black holes form is partly due to the "grid-stretching" phenomena associated with 
maximal slicing. 
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FIG. 9: Surface plots of ip at t — and t = 3 from the h adaptive simulation discussed in Sec. IV Bl where the height of each 
surface is proportional to the magnitude of (ranging approximately from 1.00 (1.00) at the outer boundary to 1.24 (3.20) at 
the origin p — z = 0a,tt = (t = 3)). Overlaid on the surfaces are the AMR grid bounding boxes — there are 4 levels of 2 : 1 
refinement at t = 0, 7 levels at t = 3, and the base (coarsest) level has a resolution of 257 x 513. 

VI. CONCLUSIONS 

In this paper we have discussed modifications of the standard Berger and Oliger adaptive mesh refinement method, 
so that the resulting algorithm can solve systems resulting from the discretization of coupled, non-linear, hyperbolic and 
elliptic equations in linear time. Moreover, as we have retained recursive time-stepping, the algorithm is still optimally 
efficient in solving systems that contain hyperbolic equations which are subject to a CFL stability condition. The 
initial application that drove development of the method was a study of critical gravitational collapse of the scalar field 
in axisymmetry |lllll2tl3^ |. This involved the approximate solution of a mixed elliptic/hyperbolic system: the coupled 
Einstein-Klein Gordon system in a certain coordinate system, with a particular choice among the overdetermined set 
of PDEs at our disposal to advance the solution in time. However, the algorithm is sufficiently general that it can 
be applied to a variety of similar systems of partial differential equations. For example, there are numerous problems 
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FIG. 10: Convergence factors 11311 for the variables ip,a,(3 p and f3 z , calculated from the fully constrained adaptive simulations 
discussed in Sec. IV Bl For second-order accurate finite differencing one would expect Q m 4. Time symmetric initial data was 
used in this case, for which the exact solution for /3 P and /3 Z at t — is 0, hence the anomalous spikes for the corresponding 
convergence factors then (i.e. Q ~ 0/0). A black hole forms at t ~ 2.5 (denoted by the dashed vertical line), after which 
significant gradients in metric functions develop due to the "grid-stretching" property of maximal slicing as the spacetime 
singularity is first approached, but ultimately avoided; this appears to be the dominant factor causing these 3 simulations (h, 2h 
and 4h) to start to depart from the convergent regime (though for any given resolution one expects departure from convergence 
to eventually occur). We cannot explain why the convergence factor is somewhat greater than four during intermediate times 
in the simulation, however this is not atypical behavior in simulations we have looked at. 



in astrophysics that need to evolve various matter equations coupled to gravity and/or the electromagnetic field, 
including cosmological structure formation, stellar evolution, supernovae, jets, accretion disks, etc. Many of these 
scenarios have a large range of spatio-temporal length scales that need to be modeled, and can benefit tremendously 
from Berger and Oliger style AMR. In some situations Newtonian gravity is sufficient to accurately describe the 
physics, and due to the linearity of the Poisson equation, modifications to Berger and Oliger as described here are 
not strictly necessary. However, some of the most interesting astrophysical events occur in regions where gravity is 
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FIG. 11: Convergence factors l|15|l for the assumed residual quantity ij) c — ipf from the adaptive simulations discussed in 
Sec. IV Bl where ip c is the conformal factor %j) from the fully constrained simulation, and ijjf is the conformal factor obtained via 
free evolution. At t = both ipf and ip c are calculated by solving the Hamiltonian constraint, hence the anomalous behavior 
in the convergence factor then. Moreover, the boundaries conditions emplo yed are only consistent with the full set of Einstein 
equations in the limit where the outer boundary position goes to infinity |32|1 . Early on in the simulation this inconsistency 
produces a difference in ip c — ipf that is of the same magnitude as the truncation error of the h simulation, and this causes 
Q^ c -ipf t° start below 4 near t = 0. However, as evolution proceeds and gravitational collapse occurs, the gradients in ip in 

the central region grow rapidly, and the truncation error component of Q^ c -^j begins to dominate. As with the results shown 
in Fig. I1UI grid-stretching effects apparently cause the decrease in Q after the black hole forms near t ~ 2.5. However, in this 
plot we can see a trend to improved convergence results at any given time when resolution is increased from Q 2h to Q . 

sufficiently strong that nonlinear effects become important, and the algorithm described in this paper could be of 
significant use in such simulations 12 . 

With regards to new applications of this algorithm in numerical relati vity , of particular interest is constrained 
evolution in 3D, which, on the basis of substantial ID and 2D evidenced Ea E3, H3, E3, H3, E3, H3, H3, H3, E3, EE 
Ml |5£| lU |62, |&|, we have long felt has great potential for the study of problems such as black hole collisions, 
critical gravitational collapse, and the structure of black hole interiors. All of the aforementioned lower dimensional 
studies made use of the symmetries in the problem, in addition to particular coordinate choices, to obtain well-posed 
coupled elliptic/hyperbolic systems, and so there is some skepticism in the community about whether constrained 
evolution can be implemented for general problems in 3D. Though recently a fully constrained 3D evolution scheme 
was proposed in |38j | , based on the Dirac gauge and spherical coordinates (the implementation presented in |38j | made 



1 




i i i i 


1 1 1 1 


— y 

— f- 


Q2h 


— 
— 

\ '*. — 










Qh 






III 


1 1 l\l 



We note that even the conformally fla t appro ximation to the field equations, which seems to be an adequate extension of Newtonian 
gravity for a certain class of problems l59ll6Cl . requires solution of a non-linear elliptic equation. 
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use of a multidomain spectral solution method). A potential disadvantage of this system is that it is not obvious 
how to generalize it beyond spherical coordinates, which are not well adapted to studying problems that are far from 
spherical symmetry. A formulation of the Einstein equations that is amenable to 3D constrained evolution and can be 
written in Cartesian coordinates was described in [64( • Both of these formulations are prime candidates for numerical 
implementation using our AMR algorithm. 
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APPENDIX A: EQUATIONS AND FINITE DIFFERENCE OPERATORS 

Here, for completeness, we list all of the equations introduced in Section [HI and the specific set of finite difference 
operators used to discretize them. 



1. Equations 



Summarizing Section [HJ the four- metric is 



ds 2 = (-a 2 + i/> 4 [((3 P ) 2 + (f3 z ) 2 ] ) dt 2 + 2t/> 4 ((3Pdp + (3 z dz) dt 

+V/ {dp 2 + dz 2 + p 2 e 2ps d<t> 2 ) 



(Al) 



The conjugate variable to the scalar field $ is II l|A8(l . and the conjugate to a is Q, All of these variables are 
functions of p, z and t. The maximal slicing condition results in the following elliptic equation for a: 
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The Hamiltonian constraint is 
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The p and z momentum constraints are 
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The definition of f2 (jBJ) gives an evolution equation for a: 



a t =2pP (pa) p2 +(3 z a, z -an 
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The evolution equation for fi is 
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The definition of II and the wave equation for <I> give 



$ t = /? p $,p + /3 2 $, z + ^n, 



(A7) 



(A8) 



and 



n f = /3Ti. p + /^n 2 + -n (a P n + 2p p p + p% 



i 



2 (pa^ 2 *,p) o2 + (aV> 2 $ . 



a 

+ ^2 



The maximal slicing condition K — gives an independent hyperbolic evolution equation for ip: 

j, = + iP, P P P + | (2f3 p p + P% + pan) . 
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2. Boundary Conditions 

The set of axis regularity conditions, applied at p = are: 
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3. Finite Difference Operators 

In this section, we write out all of the difference operators used to convert the differential equations in the previous 
section to finite difference equations. At all interior points of the mesh, the centered forms of the derivate operators are 
used, and along boundaries, backward and forward operators are used as appropriate. Kreiss-Oliger style dissipation 
is applied to evolution equations, at interior points at least two grid points inward, in the direction of the stencil, from 
any boundary. For a and £l, we linearly interpolate in p at location Ap (and optionally at 2Ap as well), using the 
values of these variables at p — and p = 2Ap (or p — 3Ap). Below, we use the notation Uij to label a point in the 
mesh corresponding to coordinate location ((i — l)Ap, z m - m + (j — l)Az (except for the coordinate variable p, where 
it is sufficient to use pi). For time derivatives, we use u™ 3 to denote the retarded time level, and u"^ 1 the advanced 
time level. All of the finite-difference operators are 2 nd order accurate. 

a. Centered Difference Operators 
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b. Forward-Difference Operators 
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c. Backward-Difference Operators 
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d. Dissipation Operators 



The following dissipation operator is applied in the p direction: 



(A25) 



and in the 2 direction: 
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where < ed < 1. 



APPENDIX B: ADDITIONAL ALGORITHM DETAILS 



This appendix contains descriptions of a few additional features of the AMR algorithm described in the paper: 
computing TE estimates using a self-shadow hierarchy, details of how the multigrid algorithm is applied to an adaptive 
hierarchy, and a description of the particular set of interpolation and restriction operators used. 



A self-shadow hierarchy is a simplification of the idea of using a shadow hierarchy to do truncation error estimation. 
A shadow hierarchy is a coarsened (usually with p s = 2 : 1 version of the main hierarchy. Both hierarchies are 
evolved simultaneously, and the function values of a given grid in the shadow hierarchy are replaced with those of 
the corresponding grid in the main hierarchy whenever the two are in sync. For example, with p t — 2 : 1, each 
time step of a shadow grid corresponds to two time steps of the main grid, and the shadow is updated every two 
main-grid time steps. A TE estimate can therefore readily be computed by comparing function values in the shadow 
with corresponding values in the main hierarchy just before the update step. 

Notice however, that within the recursive time-stepping flow of the Bcrgcr and Oliger algorithm, information for 
computing a TE estimate is "naturally" available prior to the fine-to-coarse grid injection step (see Fig. 0}. The 
coarse level £ is evolved independently of the fine level I + 1 from to to to + At c , where At c is the coarse level time 
step. Also, at i = to the level £ grid functions are restricted copies of level £+1 grid functions in the region of overlap 
O e e +1 . Therefore, prior to injection at time to + At c , the difference in an evolved variable / in levels I and £+1, within 
the region 0^ +1 , can serve as an approximation to the truncation error r s (f£+i) for / at level I + l: 13 



Therefore, for levels I > 1, one can use (|B1|) as the basis for computing truncation error estimates, without the need to 
refer to a shadow hierarchy (i.e., the main hierarchy "casts its own shadow", hence the name self-shadow hierarchy). 
This method cannot give a TE estimate for the coarsest level (1) in the hierarchy, and so we require that the coarsest 
level always be fully refined. Thus, the resolution of level 2 is chosen to match the desired coarsest resolution for a 
given problem (for the sample evolutions described in in Sec. [V] level 2 is always quoted as the base level). 

In practice, a slightly modified form of lBll is used, as we now describe. Depending upon the problem, one or more 
of tp, $, n$, Cl, and a are used in the calculation (i.e., all evolved quantities — ip is not used when the Hamiltonian 
constraint is used to solved for tp in a fully constrained evolution). Optionally, the truncation error estimate is scaled 
by the norm of the function, if the norm is larger than some constant k (k = 1 typically), and can also be multiplied 
fp , for some integer p chosen heuristically to either enhance or reduce the near- axis refinement. The TE estimate for 
a function ft at level £ is thus defined to be 



In fact, if there are only two levels in the hierarchy, then this estimate is exactly the truncation error estimate one would obtain with 



1. A Self-shadow Hierarchy for Computing Truncation Error Estimates 
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where it is implied that ft is only defined in the overlap between levels £ and I — 1, fg is restricted to the resolution 
of level ^ — 1 prior to subtraction, and the result is then interpolated back to the resolution of level £. Typically, we 
use p = 2 for Cl, p = 1 for a, and p = for the other variables. The TE estimate for the level is defined to be 



where the sum is taken over the desired subset of variables listed above. 

Optionally, the TE estimate calculated in i|B3|l is further smoothed (using simple averaging over a 5-by-5 square 
cell of points), and/or scaled by a quantity Hi > 1 in the region of overlap between levels £ and £ + 1. Hi therefore 
provides a kind of "hysteresis" to the truncation error estimation process: when the TE estimate in a region of level 
£ grows above r max , that region is refined; however, for the region to be unrefined at a later time, the TE estimate 
needs to drop below T max /i?£ there. In most of the simulations we keep Hi = 1, though occasionally it has proven 
useful to set it to around 5 — 10. 



The FAS multigrid algorithm, with ^-cycling, that we use to solve elliptic equations on a grid hierarchy such as 
that shown in Fig. [21 is based on Brandt's multi level adaptive (MLAT) scheme ^6], and is similar to that used in 
|25|. To simplify the algorithm, we require that p s = 2 q for some integer q; then the AMR hierarchy can easily be 
extended to incorporate the multigrid levels, which have a refinement ratio p mg of 2:1. When building the multigrid 
hierarchy, each AMR grid is individually coarsened by factors p mg (i.e. factors of 2) until either a) one dimension of 
the grid is smaller that the minimum allowed, or b) the coarsened grid can be "absorbed" into a larger grid at that 
level in hierarchy. With this type of hierarchy the y-cycle begins with relaxation of the finest level grids only. Then 
as the finer levels are coarsened they are absorbed, if possible, into coarser levels. This method of relaxation on an 
adaptive hierarchy is in contrast to the Fast Adaptive Composite Grid Method 65] , though we do not believe that 
the particular details of how the elliptic equations are solved have much bearing on the algorithm described here for 
dealing with coupled elliptic/hyperbolic equations. 

A couple of minor points regarding the details of the multigrid solution are worth mentioning. First, inter- level 
operations, such as restriction, computing coarse-grid corrections, etc., are only performed in the region of overlap 
between the two levels (which is always the region of the fine level, given the kind of hierarchies that are produced 
by B&O AMR). Second, the manner in which the relaxation sweep proceeds over a level is modified, to account for 
possible grid overlap 14 , as follows. During the relaxation sweep, the variables at a given coordinate location and 
level are relaxed only once, regardless of the number of grids encompassing that location. This is crucial in order to 
preserve the smoothing properties of the relaxation scheme. We use a mask function to enforce this requirement of 
a single update per grid point. The mask is initialized to zero on all grids at that level, prior to a relaxation sweep. 
Then, on a given level, a sweep is applied, in turn, to each grid in the level, but only variables at points where the 
mask is equal to zero are modified. After the sweep is complete on a given grid, the mask is set to one throughout the 
grid, and the mask and other grid functions are copied to overlapping grids at the same level. Therefore, subsequent 
relaxation sweeps on adjacent grids skip over points that have already been relaxed. This communication step, in 
addition to enforcing a single update per point, ensures that grid functions are numerically unique at all physical 
grid locations 15 , which is important for preserving the convergence properties of multigrid. Also, although Dirichlct 
boundary conditions are used at interior boundaries (i.e. those not abutting boundaries of the computational domain) 
of any single grid, the communication ensures that points interior to a union of grids are ultimately updated using 
the PDEs, even if they lie on the boundary of some grid in the overlap region. 

With regards to the performance of this multigrid scheme on a general adaptive hierarchy, there are two situations 
of relevance where performance could suffer, compared to the single grid multigrid algorithm. The first occurs when, 
at some level down (coarser) in the multigrid hierarchy, one or more grids in a connected union of grids is a "coarsest 
grid" , and hence needs to be solved "exactly" 16 — see Fig. Experimentation showed that the entire union needs 



We require that grids at the same level must overlap when spanning a connected region of the computational domain. In other words, it 
is not sufficient for grids to merely "touch" at a common boundary between them. This requirement simplifies the relaxation subroutine 
so that it can operate locally on a grid-by-grid basis, without needing to communicate adjacent information. However, as described in 
the text, the communication step then needs to be shifted to other parts of the algorithm. 

Note that such a communication step is also performed after relaxation of evolved variables, during the Crank-Nicholson iteration. 
Here, "exactly" means that the residual on the coarsest grid is reduced by several orders of magnitude by relaxation (it is usually not 
necessary to solve the coarse-grid problem to within machine precision). 




(B3) 



2. Multigrid on an Adaptive Hierarchy 
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FIG. 12: A hypothetical multigrid grid configuration that will adversely affect execution speed. In the figure we depict three 
grids, gl,g2 and g3, several levels down in the grid hierarchy (i.e., after several coarsening steps have already been performed). 
Grid g3 cannot be coarsened any further, while grids gl and g2 can, and ideally should be coarsened further, to maintain the 
speed of the algorithm. However, because g3 overlaps the other two grids, this entire level must be considered a coarsest level, 
and solved "exactly". 

to be solved exactly in that situation; i.e. it is not sufficient to solve the equations exactly on the coarsest grids, then 
proceed down the y-cycle on the remaining grids. If the union of grids consists of a relatively small number of grid 
points, then such a situation will not be a problem; otherwise, there will be a significant slow-down of the code, for 
the speed of an exact solve suffers dramatically as the number of unknowns increase. To date, we have been able to 
avoid this potential speed bottleneck by using a more simplistic clustering method that does not produce grid-overlap, 
as discussed in the following section. 

The second situation where performance suffers is when the TE estimate requires long, skinny rectangular regions 
to be refined. This does occur with the more prolate initial data configurations that have been studied 0]. What 
happens then, is that such an elongated grid can not be coarsened very much along the larger grid dimension before 
the smaller dimension has reached the smallest allowed size. Again, this results in a relatively large number of points 
on the coarsest grid where the solution needs to be obtained exactly. As of yet this problem has not been addressed. 

3. Clustering 

We have incorporated two clustering routines into the code. The first, written by Reid Guenther, Mijan Huq and 
Dale Choi [(5(| , is based upon the signature- line method of Berger and Rigoutsos [63 ■ The second is a simple routine 
that produces single, isolated clusters — each isolated region of high TE is surrounded by a single cluster, and then all 
clusters within a certain distance of each other are merged together into an encompassing cluster. For the problems 
studied so far, the isolated cluster method turns out to be almost as efficient as the signature-line method. Therefore, 
since efficiency is not an issue, the isolated cluster method is preferable, because it avoids one of the potential speed- 
bottlenecks of the multigrid algorithm discussed in Section IB~2"1 furthermore, as mentioned in Section IB~51 minimizing 
cluster overlap helps reduce high-frequency noise problems. 

A clustering issue that needs to be dealt with in this code is that the resultant grid hierarchy must be compatible 
with the multigrid solver. This places two restrictions on cluster sizes and positions. First, an individual grid must 
have dimensions that can be factored into x m i n 2™, where x m i n is one of the smallest, allowed grid dimensions, and 
n is a non-negative integer. Second, if several grids overlap, then their relative positions must be such that the 
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common grid points align on all possible levels of a multigrid hierarchy. Specifically, if a union of overlapping grids 
can collectively be coarsened m times in the multigrid hierarchy, then the relative offsets of grid origins on the finest 
level must be multiples of 2 m grid points. These requirements are enforced after the initial clustering algorithm is 
called, by modifying the returned cluster list accordingly. This gives more flexibility to experiment with different 
clustering routines, which consequently do not need to be aware of the alignment issues. 

To conclude this section we mentioned a couple of additional options that have been implemented in the post- 
clustering routine. They are adding "ghost zones" between adjacent, touching clusters, so that both the multigrid 
and evolution relaxation sweeps correctly solve the system of equations in a domain given by the union of grids 
at a given level; and optionally moving or shrinking clusters, if necessary, to prevent them from touching parent 
boundaries 17 , which helps to avoid instabilities that occasionally occur in such situations. 

4. Interpolation and Restriction Operators 

Here we state the restriction and interpolation operators used in the AMR code. Simple injection is used to restrict 
a fine grid to a coarse grid during the injection phase of the AMR algorithm, and when computing the TE estimate. 
A fourth order (bi-cubic) interpolation scheme is used to initialize newly refined fine grids (or regions thereof) from 
the encompassing coarser grid. The scheme proceeds by first interpolating every row of the coarse grid to the fine 
grid (i.e. every p s th row of the fine grid is filled in), then all the remaining points on the fine grid are computed via 
interpolation, column-by-column. The multigrid routine uses half-weight restriction when transferring from a fine to 
coarse grid, and linear interpolation for the coarse to fine transfer. 

5. Initializing the Grid Hierarchy 
Fig. EH below contains a pseudo-code description of the steps used to initialize the grid hierarchy. 

6. Controlling High-Frequency Grid-Boundary Noise 

An issue that needs to be dealt with in a Berger & Oliger style AMR scheme is controlling high-frequency solution 
components ("noise") that may occur at parent-child grid boundaries. For a second order accurate finite-difference 
scheme, the second derivatives of grid functions are typically not continuous across the boundaries after child to 
parent injection. This potential source of high-frequency noise on the parent level is rather efficiently eliminated by 
the Kreiss-Oliger dissipation filters we incorporate into our finite differenced evolution equations. 

In certain situations we have found that high-frequency noise also develops on child grids, within a grid point or 
two of the AMR boundary (in particular near the corners of the grid, or places where two grids overlap). This noise 
is not as easily dealt with, as the Kreiss-Oliger filter acting normal to the boundary is only applied a distance three 
points and further away from the boundary. The source of this noise appears to be the parent-child interpolation 
scheme used to set the boundary values, and in general the interpolation method must be tailored to each variable in 
order to reduce the noise to an acceptable level. For out current model, we use the following interpolation method. 
For all evolved variables (a, f2, $, 11$ and ip) we use linear interpolation in time from the parent level to set boundary 
values on the child grid at points coincident with parent grid points. This is followed by fourth-order interpolation 
in space (as described in Section TB 4|) for the remaining boundary points. Furthermore (see Fig. 1141 and IT5f >. after 
each step of the Crank-Nicholson iteration we reset f2 and a in a zone two grid points in from AMR boundaries with 
values obtained either 1) by fourth order interpolation using function values from the boundary and three additional 
points inward from this zone, or 2) via bilinear interpolation at "corner" points, i.e. those points that are a single 
cell width away from two boundaries. This technique for a and f2 was discovered after quite a bit of experimentation 
with different interpolation schemes, and is quite effective in reducing the level of noise at the grid boundaries. 

For the elliptic variables (a,(3 p , f} z and optionally ip), prior to a Crank-Nicholson evolution cycle, we reset these 
variables on AMR boundaries at points unique to the grid (in between points coincident with parent level points — 



In principle this should never occur if one adds a buffer zone about the region of high truncation error. However, because of the grid 
shuffling performed to obtain a hierarchy acceptable to multigrid, a grid could be extended to touch a parent boundary. With the option 
enabled to prevent this, the grid will be reduced rather than extended to fit into the multigrid scheme. This comes at the expense of 
not obtaining "optimal" zones about the region of high truncation error. 
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t:=0; 

initialize the grid hierarchy with 2 levels, each covering the entire domain; (so if = 2) 
repeat 

call set_initial_data() ; (see below) 

call set_past_t_data_lst_order () ; (see below) 

call single_step(l) ; (see Fig!^ 

regrid the entire hierarchy using truncation error estimates 
computed in the previous step; (thus possibly changing if) 
reset t to t := while retaining the current hierarchy structure; 
until if — i P f or maximum number of refinement levels reached 

call set_initial_data() ; 

call set_past_t_data_lst_order () ; 

call single_step(l) ; (evolve hierarchy forwards in time one coarse step) 
call flip_dt; (see below) 

call single_step(l) ; (evolve hierarchy backwards in time one coarse step) 
call flip_dt; 

call set_initial_data() ; 

(done computing initial data and hierarchy) 

subroutine set_initial_data() 

initialize hyperbolic variables over levels [1..^/] with freely-specifiable data; 

solve the elliptic equations over levels [l..£/]; 
end of subroutine set_initial_data 

subroutine flip_dt() 

for each elliptic variable fi(t): fi(t + Ati) := 2/ x (t) - fi(t - Ati) ; 
do I := 2 to if 

for each elliptic variable fi(t): f t (t + p t At t ) := 2f e (t) - f e (t - p t At t ) ; 
At e := -At/; 
end do 

end of subroutine flip_dt 

subroutine set_past_t_data_lst_order () 

do £ := 1 to if 

for each elliptic variable fi(t): fi(t — p t Ati) :— fe{t); 
end do 

end of subroutine set_past_t_data_lst_order 

FIG. 13: A pseudo-code description of the steps we use to initialize the grid hierarchy. The repeat loop is used to calculate the 
hierarchy structure at t = 0. Then, to initialize past time level data for elliptic variables using linear extrapolation (which is 
done in flip_dt()), the entire hierarchy is evolved forwards, then backwards in time by single coarse level time steps (alternatively 
we could evolve backwards, then forwards in time here — the results would essentially be the same). The idea behind this last 
step is that since the evolution of the hyperbolic variables is driving any change in the elliptic variables, we can use the results 
of a small evolution step to provide a better estimate of past time level information than first-order extrapolation of the solution 
at t = 0. In principle, this step can be iterated if need be, though we found that a single step is sufficient (though not always 
necessary depending on the free initial data) to obtain close to second order convergence of the final solution. 
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Interpolation Method 
□ (1) linear in time from parent grids 



■ (2) A th order in space using points from (1) 



• (3) 4 th order in space using a point from (1) or (2) 
and 3 interior points normal to the boundary 

X (4) as (3), though utilizing points in (3) instead of 

interior points 
o (5) bilinear in space using points from (1) and (4) 



FIG. 14: An illustration of the interpolation method used for a and Q during AMR evolution. In the figure we assume that 
the spatial refinement ratio is 2 : 1, and that all four grid boundaries are AMR boundaries. Points labeled by (1) and (2) are 
set once prior to the Crank-Nicholson (CN) iteration, while points labeled by (3), (4) and (5) are reset after every CN step (see 
the pseudo-code in Fig. I15H . Points not explicitly labeled are "interior" points, and are evolved. 

subroutine interp_interior_AMR_bnd(grid function f [1 . . Nrho , 1 . .Nz] ) 
for i=3 to Nrho-2 do 

f (i,2)=0.4*f (i,l)+2.0*f (i,4)-2.0*f (i,5)+0.6*f (i,6) 

f (i,3)=0.1*f (i,l)+2.0*f (i,4)-1.5*f (i,5)+0.4*f (i,6) 

f (i,Nz-l)=0.4*f (i,Nz)+2.0*f (i,Nz-3)-2.0*f (i,Nz-4)+0.6*f (i,Nz-5) 

f (i,Nz-2)=0. l*f (i,Nz)+2.0*f (i,Nz-3)-1.5*f (i,Nz-4)+0.4*f (i,Nz-5) 
end do 

for j=3 to Nz-2 do 

f (2,j)=0.4*f (l,j)+2.0*f (4, j)-2.0*f (5,j)+0.6*f (6, j) 
f (3,j)=0.1*f (l,j)+2.0*f (4,j)-1.5*f (5,j)+0.4*f (6,j) 

f (Nrho-1 , j ) =0 . 4*f (Nrho , j ) +2 . 0*f (Nrho-3 , j ) -2 . 0*f (Nrho-4 , j ) +0 . 6*f (Nrho-5 , j ) 
f (Nrho-2 , j ) =0 . l*f (Nrho , j ) +2 . 0*f (Nrho-3 , j ) -1 . 5*f (Nrho-4 , j ) +0 . 4*f (Nrho-5 , j ) 
end do 

f (2,2)=(f (1,1) +f (3,3)+f (l,3)+f (3,l))/4 

f (Nrho-1 , 2) = (f (Nrho , 1 ) +f (Nrho , 3) +f (Nrho-2 , 1) +f (Nrho-2 , 3) ) /4 
f (2,Nz-l)=(f (l,Nz)+f (3,Nz)+f (l,Nz-2)+f (3,Nz-2))/4 

f (Nrho-1 , Nz-1) = (f (Nrho , Nz) +f (Nrho , Nz-2) +f (Nrho-2 , Nz) +f (Nrho-2 , Nz-2) ) /4 
end of subroutine interp_interior_AMR_bnd 

FIG. 15: A pseudo-code description of part of the interpolation method used for a and Cl during AMR evolution. For simplicity, 
in this subroutine we assume that all four boundaries are interior to the computational domain boundaries. The set of points 
altered here correspond to (3), (4) and (5) in Fig. 1141 and the interpolation operators used are independent of the spatial 
refinement ratio (as opposed to points (1) and (2) in Fig. 1141 . 

those labeled by (2) in Fig. I14f> using fourth order interpolation from the remaining points on the boundary. This 
overwrites the values set by coarse-grid corrections (CGCs) during the most recent multigrid solve that involved 
coarser levels. The reason for doing this is as follows. In multigrid, CGCs typically introduce high-frequency noise on 
the finer level, while the subsequent post-CGC relaxation sweeps smooth out this noise. However, since no relaxation 
is applied on AMR boundaries, some form of explicit smoothing is required — the above described fourth order 
interpolation provides this smoothing mechanism. 

Another stage in the algorithm where high-frequency noise can creep into the solution is during the regridding 
phase, if the refined region on a given level expands. Then, within the part of a new grid overlapping the old refined 
region, grid functions will be initialized by copying data from an old fine grid, while on the remaining part of the new 
grid the grid functions will be initialized via interpolation from a parent grid. Sometimes, at the interface between 
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the copied/interpolated data, tiny discontinuities are introduced. The grid functions are then easily smoothed by 
applying a Kreiss-Oliger filter to them (at all time levels involved) after regridding. 



